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Abstract 

Using equilibrium and nonequilibrium molecular dynamics simu- 
lations, we determine the Kapitza resistance (or thermal contact re- 
sistance) at a model liquid solid interface. The Kapitza resistance (or 
the associated Kapitza length) can reach appreciable values when the 
liquid does not wet the solid. The analogy with the hydrodynamic 
slip length is discussed. 

1 Introduction: bulk and interfacial trans- 
port coefficients 

One of the important issues of nonequilibrium statistical physics is the deter- 
mination of transport coefficients, which relate the flux of conserved quanti- 
ties to the gradient of the corresponding affinities. One of the simplest ex- 
amples is the thermal conductivity defined macroscopically (in the absence 
of concentration or velocity gradients) by 

je = -AVT = LeV^. (1) 
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Here Lg is the Onsager coefficient associated with energy transport, 1/T = 
dS/dE is the affinity relative to energy. 

Molecular dynamics, using either equilibrium or nonequilibrium methods, 
has proven a powerful tool for determining such coefficients, and, in fact, the 
only practical one when dealing with dense liquids. Equilibrium methods 
for systems with continuous potentials, pioneered by Levesque and cowork- 
ers are based on the use of Green-Kubo formulae p|. For the thermal 
conductivity, the formula reads: 



where je is a microscopic energy current that can be computed directly from 
the particle velocities and positions. The brackets refer here to an equilibrium 
average on the microscopic fluctuations at equilibrium. Green-Kubo calcu- 
lations have the advantage that an equilibrium simulation can be used to 
determine simultaneously several transport coefficients. On the other hand, 
rather long simulations are required to obtain converged results for the time 
correlation function, and the determination of the 'plateau' of the Green- 
Kubo integral (as a function of the upper integration limit) is often difficult, 
since at long times the noise in the correlation function tends to produce a 
large numerical uncertainty. 

A useful alternative for the determination of transport coefficients is pro- 
vided by nonequilibrium molecular dynamics (NEMD) methods In 
such methods, the system is driven into an equilibrium state using a mod- 
ified MD algorithm, and the currents and gradients are measured directly 
in this state. Formally, these methods can be shown to be equivalent to the 
Green-Kubo equilibrium method as long as one remains in the linear response 
regime. 

In both types of approach, one usually is interested in bulk transport coef- 
ficients, which are properties of the infinite, homogeneous system. Much less 
attention has been devoted to transport processes across interfaces between 
two bulk phases. At the macroscopic scales, the flux of a conserved quantity 
must be conserved across the interface. Continuum theories, however, usually 
assume continuity not only of the fluxes, but also of the affinities. The latter 
assumption, however, has no theoretical justification, and is only known to 
work 'in practice' at the macroscopic scale. A more consistent description 
of transport across an interface, in the spirit of usual hydrodynamic descrip- 
tions, allows for a jump in affinities between the two phases, proportional to 




(2) 
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the flux of tlie conserved quantity. In tlie case of energy transport, such a 
phenomenological description can be translated mathematically as 

ji.-ni2 = :^(T(^)-T(^^ (3) 
±1k 

To conform with usual notations, we have used the jump in temperature 
rather than the jump in affinity 1/T. Equation |^ defines the Kapitza resis- 
tance Rk (or thermal contact resistance) of the interface, whose existence 
was pointed out by Kapitza in the context of liquid Helium physics in |^ . A 
more transparent definition can be obtained by defining the Kapitza length 

Ik = RkX (4) 

where A is the thermal conductivity of one of the phases (say phase (1)). 
Ik is simply the thickness of material (1) equivalent to the interface from a 
thermal point of view (see figure 0). The Kapitza resistance has been mostly 
studied in the case of superfiuid Helium in contact with a solid, or in the case 
of grain boundaries between pure crystals 0. In both cases A is large, so 
that the Kapitza length is large and measurable effects can be expected. Size 
effects on the thermal conductivity of composites are also an indirect way of 
assessing the thermal properties of interfaces. In the case of normal liquids 
in contact with a solid, however, it can be expected, on purely dimensional 
grounds, that Ik should be of the order of a few molecular sizes, and therefore 
not measurable. 

In this paper, we show that Ik is indeed very small in the -most usual- 
case of a wetting (contact angle smaller than 7r/2) liquid-solid interface. In 
the nonwetting case, however, the coupling between the solid and the liquid 
is strongly reduced, and Ik can be increased considerably - reaching up to 50 
molecular sizes. 

Before we discuss in more detail the determination of Ik for a model 
liquid-solid interface, it is interesting to point out an analogy with another 
interfacial transport coefficient that has attracted considerable attention in 
the recent years, namely the 'slip length'. The definition of the slip length 
completely parallels that of the Kapitza length. In the bulk, the phenomeno- 
logical equation relates the momentum fiux (stress) to the velocity gradient 
(up to a factor of density, the velocity is the affinity relative to momentum). 
The transport coefficient of interest is the shear viscosity t], which in a simple, 
Couette fiow situation is defined by 

CTxz = (5) 
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Figure 1: Schematic illustration of the Kapitza length (left) and of the slip 
length (right) at a liquid solid interface. 
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where x is the flow direction and z the direction of the velocity gradient 
(see flgure |l]). At the interface, the tangential velocity is not continuous in 
general. Since Vx is zero in a solid at rest, the equivalent of equation ^ reads 

(Jxz{z = 0) = nvxiz = 0). (6) 

Combining the two equations | and one obtains a boundary condition for 
the velocity fleld 

^ = -Vx{z = 0) = -Vx{z = 0). (7) 

OZ 1] 

Equation |^ deflnes the slip length 6, which is the distance over which the 
velocity fleld must me extrapolated inside the solid to obtain a zero velocity 
(see flgure 0). In macroscopic hydrodynamics, 6 is usually assumed to van- 
ish, which amounts to a 'stick' boundary condition at the interface. However, 
with the development of accurate surface force measurement at the nanome- 
ter scale, it has become possible to measure nonzero values of 6. From a 
theoretical standpoint, it has been possible to determine S (or rather the 
coefficient k in equation ^) using a Kubo like formula: 

where F^^^ is the instantaneous force along x exerted by the liquid on 
the solid across the interface. It has also been shown, both theoretically 
and experimentally, that the wetting properties are an important ingredient 
in determining 6. Everything (lattice constants, densities) being otherwise 
equal, changing the solid-hquid interaction (and therefore the wettabihty) 
changes 6 by orders of magnitude. Large (50 nm or more) values of 6 have 
been measured for water on hydrophobic surfaces p|, |^, [10|, |Tl| . 

Our model of the liquid-solid interface will be identical to the one used 
for our previous study of slip lengths. The interactions between species (i) 
and (j ) are of the Lennard- Jones form 



Vij{r) = 4e 



(9) 



in which, for convenience, the strength of the attractive term is modulated 
by a coefficient Qj, keeping e and a identical (obviously the same effect could 
be achieved by introducing different values of e and a; in this form, emphasis 
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Figure 2: Transverse cut of a typical simulation cell. In order to avoid any 
possible problem with dissolution of solid atoms, those atoms are loosely 
bound to their equilibrium position through weak harmonic springs. 

is put on the attractive part which is essential in controlling the wetting 
properties). For the liquid-liquid interactions, we set cn = 1.2. For the FCC 
solid, C22 = 1, and the sohd atoms are tethered to fixed lattice sites at a 
constant density, PsCr^ = 0.9, using weak harmonic springs. All interactions 
are cutoff and shifted at 2.5cr. The interfacial properties of the (100 FCC 
solid with the liquid have been studied in detail in reference |T^. By varying 
the coefficient C12, it is possible to vary the contact angle 6, up to 6* = 140° 
for C12 = 0.5. 

In the next two sections, we describe the methods used for determining the 
Kapitza resistance of this model interface. Results are presented in section 

i 

2 Nonequilibrium molecular dynamics simu- 
lations 

The setup used for computing the Kapitza resistance is illustrated in figure 
^ A slab of liquid (typically of thickness between 30cr and 40ct) , parallel to 
xOy is sandwiched between two layers of FCC solid of thickness 10 — 20cr. 
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Figure 3: Temperature profile obtained for Ci2 = 0.5 

The lateral dimensions of the cell, depending on the simulation, have been 
varied between 10 and 20a. Periodic boundary conditions are applied in the 
X and y directions. Liquid densities are typically of order pia^ = 0.8 — 1, 
and the reference temperature is always T = le/ks- The normal pressure P 
is obtained from the average force exerted by the liquid on the solid; all the 
simulations have been carried out at low values of the pressure, \P\ < O.le/ cr^. 

A temperature gradient is generated throughout the cell by coupling the 
upper and lower atomic layers of the solid to two thermostats at different 
temperatures. The thermostats use simple rescaling of the velocities for 
the group of atoms under consideration. When the two temperatures are 
different, a net energy flux Je in the z direction results, which is computed 
by averaging the kinetic energy added (or removed) by each thermostat per 
unit time and surface. Temperature proflles are obtained from the local 
kinetic energy density. 

A typical temperature proflle obtained under these nonequilibrium con- 
ditions is shown in flgure In the 'bulk' liquid and solid, the temperature 
varies essentially linearly with position, in agreement with the predictions 
of Fourier's law. Note that a slight curvature is perceptible in the liquid 
phase, which reflects the dependence of A on temperature. From the slope of 
the linear parts, the following values of the thermal conductivity in the bulk 
phases are obtained (all results are given in Lennard- Jones units, e, a, and 
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Xi{T = 1,P = 0) = 7.8 ±0.5 ■ X,{T = l,p = 0.9) = 8 ± 1 (10) 

The liquid has a density pscr^ = 0.88 at T = 1 and P = 0. The uncertainty 
on the resuhs has two origins: first, the energy current is the nonzero average 
value of a strongly fluctuating quantity (i.e. the power extracted by the ther- 
mostats), and is therefore difficult to determine accurately. Our simulations, 
with a time step h = 0.01 and typically 10^ steps, yield this quantity with 
an accuracy of about 5%. In the solid, the limited thickness makes the de- 
termination of the temperature gradient inaccurate, so that the uncertainty 
is somewhat larger than in the liquid. 

Another issue concerning the solid thermal conductivity is the comparison 
of the phonon mean free path iph with the thickness of the solid slab. Using 
As ~ iphUCy/3 where u is the sound velocity (obtained from the equation 
of state) and the heat capacity per unit volume, one finds iph — 2.5a. 
Therefore the thickness of the solid slabs is of order 3 — Mph, which means 
that ballistic effects in the transport of heat should be negligible. Indeed, we 
confirmed by some simulations using thicker solid boundaries {20a) that the 
results were not sensitively affected by this parameter. 

From the profile shown in figure |^ the Kapitza resistance is easily obtained 
by dividing the temperature jump (computed using the linear fits to the 
density profiles) by the energy fiux. Note that a slight dependence of Rk on 
T is visible in figure |^. with a larger temperature jump on the 'warm' side. 
All results discussed in section ^ will correspond to the average values of Rk 
between the 'warm' and 'cold' sides, and therefore to a nominal temperature 
T=l. 



3 Equilibrium determination of Rx 

A completely different route to the Kapitza resistance is provided by a Green- 
Kubo formula similar to that used for the slip length, equation ||. This 
formula, which to our knowledge was first written by Puech et al 0, reads 

Here q is the energy flux across the interface of surface S, which is macroscop- 
ically defined by the surface integral q(t) = / / j£;(t) ■ dS. Microscopically, 
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q{t) is easily computed from the work per unit time exerted by the atoms 
belonging to the solid phase on those belonging to the liquid phase, namely: 



<lit)= E E F^.-v. (12) 

illiquid jgsolid 



Equation |TT] can be justified through the following heuristic derivation. 
Let us consider a solid and a liquid coupled together through an interface of 
surface S. The total system (solid+liquid) is isolated, at equilibrium. For 
convenience, we will assume that the solid constitutes a thermostat (infinite 
heat capacity) at fixed temperature Tq, while the heat capacity of the liquid 
part is denoted by Cy- If Ei{t) is the internal energy of the liquid phase at 
time t, and if one assumes a Langevin evolution of this variable, one may 
write 

^ = gW = ;^(ro-TKt)) + /(t) (13) 

where Ti refer to the instantaneous temperature of the liquid. As usual in 
a Langevin description, we have separated the evolution of the slow vari- 
able (the energy) into a systematic contribution which would correspond 
to a macroscopic evolution and a random, 6 correlated part f{t). Writing 
(fit) fit')) = ^5it - t'), X = Eiit) - {El) = CviTiit) - To) one may rewrite 
|13| in the standard Langevin form 

^ = -aXit) + fit) (14) 



with a = n^Q^ ■ Equation is a standard Langevin equation for the variable 



X. The usual relations follow for the correlation functions 

(X(t)X(0)) = ^exp(-a|t|) (15) 

mm = -^(^W^(0)) = - ^exp(-a|t|) (16) 

In the thermodynamic limit (Cy oo) the second term vanishes, so that 
one may write: 

7 = 2/ (g(t)g(0))cit (17) 







Using the relations (X^) = y' which follows from equation p!5| at t = 



and the equilibrium distribution for X which yields the standard fiuctuation 
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Figure 4: Time correlation functions for q[t) for Cu 
(the surface of the cell is 5* = 390cr^). 



0.5, ci2 = 0.8, ci2 = 1 



formula (X^) = {6E'^) = ksT^Cy, the Kubo relation ^ for Rk is obtained 
immediately from 

As usual [0, this line of reasoning can be extended to the 'generalized 
Langevin' case where the coefficient a is replaced by a memory kernel, and 
f{t) is not a white noise. The Kubo formula is not modified, and the macro- 
scopic transport coefficient being related to the time integral of the memory 
Kernel. The Langevin approach just consists in assuming a separation of 
time scales, so that the kernel can be replaced by a 5 function with the ap- 
propriate intensity. Note, however, that the Kubo formula ^ holds only 
after the thermodynamic limit Cy — > oo has been taken. Hence, for a fi- 
nite system, the running integral 2 Jq {q{t)q{0))dt will consist of two parts: 
a first part yielding to a plateau value which contains the information on 
the Kapitza resistance, and a long time decay to zero. Typical correlation 



functions obtained in MD simulations (using the microscopic definition |T2|), 
and the corresponding time integrals, are shown in figures ^ and The 
plateau of the time integral can be used for a determination of the Kapitza 
resistance. 
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Figure 5: Running integral of the time correlation function {q{t)q{0)) ob- 
tained for Ci2 = 0.5 and P = 0.008. A value of Rk = 6.9 results, to be 
compared with the value Rk = 6.7 obtained from a nonequilibrium simula- 
tion. 

4 Dependance of Rk on wetting properties 

It has been shown previously (see |]12|) that the slip length at the liquid solid 
interface was dependent on the wetting properties. The ability of the solid 
to transfer momentum to the liquid is weaker when the liquid does not wet 
the solid, since in this case the liquid density is depleted in the vicinity of 
the solid wall. Obviously the same kind of effect could be expected in the 
Kapitza resistance, with a smaller rate of energy transfer when the liquid- 
solid interaction is weak. In order to check this expectation, we have studied, 
using both the nonequilibrium and the equilibrium approach, the Kapitza 
resistance of our model interface as a function of the interaction coefficient 
Ci2, which also determines the contact angle. All simulations are performed 
at a pressure close to P = 0. The results are shown in figure ^, where it can 
be seen that the Kapitza length is indeed a decreasing function of Ci2, with 
values close to a molecular size in the wetting case (ci2 = 1), and values up 
to about 50 molecular sizes in the nonwetting case (ci2 = 0.5). As a rule of 
thumb, Ik appears to be of the same order of magnitude as the slip length 



obtained in reference 12 
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Figure 6: Kapitza length as a function of the interaction parameter Ci2. All 
simulations are performed in the range of normal pressures 0.2 > P > —0.2, 
where we have checked that the pressure dependence of is weak. The error 
bars reflect the combination of several independent simulations (equilibrium 
and nonequilibrium) for determining each point. The heat conductivity of 
the liquid has been used to define the Ik from the Kapitza resistance. 
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5 Summary 



In this paper, we have discussed the principles underlying the definition of 
interfacial transport coefficients, and shown, in the case of the Kapitza resis- 
tance, how such coefficients can be extracted from equilibrium or nonequi- 
librium simulations. Our results indicate that relatively large values of the 
Kapitza length can be obtained using standard (as opposed to superfiuid) 
liquids, provided the liquid is not wetting the solid surface. Experimentally, 
it seems very difficult to measure directly this transport coefficient, as has 
been done for the slip length. However, Rk could probably be investigated 
by monitoring the heat conductivity of a porous matrix filled with liquid. 
Interfacial properties can then be tailored by treating the internal surface 
of the pores (e.g. grafting of short hydrocarbon chains to make the surface 
hydrophobic [Q)- 
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